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Abstract. 

Smoothing over structures in general relativity leads to a renormalisation of 
the background, and potentially many other effects which are poorly understood. 
Observables such as the distance-redshift relation when averaged on the sky do not 
necessarily yield the same smooth model which arises when performing spatial averages. 
These issues are thought to be of technical interest only in the standard model of 
cosmology, giving only tiny corrections. However, when we try to calculate observable 
quantities such as the all-sky average of the distance-redshift relation, we find that 
perturbation theory delivers divergent answers in the UV and corrections to the 
background of order unity. There are further problems. Second-order perturbations 
are the same size as first-order, and fourth-order at least the same as second, and 
possibly much larger, owing to the divergences. Much hinges on a coincidental balance 
of 2 numbers: the primordial power, and the ratio between the comoving Hubble 
scales at matter-radiation equality and today. Consequently, it is far from obvious that 
backreaction is irrelevant even in the concordance model, however natural it intuitively 
seems. 
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1. Introduction 

The standard Friedmann-Lemaitre- Robertson- Walker (FLRW) models of cosmology can 
account for all observations to date with just a handful of parameters - a resounding 
success. These are the simplest reasonably realistic universe models possible within 
general relativity: homogeneous, isotropic, and flat to a first approximation, with a scale- 
invariant spectrum of Gaussian perturbations from inflation added on top to account 
for structures down to the scale of clusters of galaxies. Were there a physical motivation 
for the value of the cosmological constant and the associated coincidence problem we 
would have a nearly complete understanding at nearly all epochs. 

The coincidence problem has led some to speculate that the evidence for a 
cosmological constant might in fact be a misunderstanding of the standard model, and 
that complex relativistic 'backreaction' effects are present when structure forms. This 
might be important in a variety of ways [1-9] . The canonical mechanism might arise by 
altering the dynamics of spacetime when we smooth over structure to form a spatially 
homogeneous and isotropic 'background' [10-18]. In the standard picture, this would be 
represented by a mismatch between the background model at early versus late times - 
in the standard model they are assumed to be the same. Alternatively, smoothing 
observations over the sky might lead to a significantly different cosmology than a 
corresponding spatial smoothing process (sometimes known as the fitting problem [19]). 
Some authors have argued that it is specifically an infra-red divergence which appears 
at second-order that gives the possibility of the backreaction effect which can mimic 
dark energy [6, 7, 11]. 

That such mechanisms might lead to significant effects - let alone be the cause of 
apparent acceleration - has been quickly dismissed by many in the community [20-36] . 
We can understand the argument why backreaction must be small as follows. 

The fields that are propagating on the background must alter its dynamics through 
the non-linearity of the field equations. According to observers at rest with the 
gravitational field, the potential itself is small everywhere outside objects less dense 
than neutron stars, and so we write g a b = gfb + 9db where g^f ~ &9if where 
$ ~ 10~ 5 <C 1. The affine connection or Ricci rotation coefficients determine the 
dynamics of the spacetime, and are generically C?(<9$); the field equations are 0(<9 2( 3>). 
A perturbation of wavelength A = a/k, where k is the comoving wavenumber and a 
is the scale factor, much less than the Hubble scale can give rise to large fluctuations 
in the field equations 0(k 2 &), even though the change from the background metric is 
small. Such terms describe density fluctuations, and these can be large even though 
the metric potentials are small. After averaging over a large domain, however, all first- 
order quantities vanish by assumed homogeneity and isotropy of the perturbations. This 
picture predicts the CMB temperature anisotropics and the matter power spectrum to 
all required precision. At second-order the Einstein tensor can only contain terms as 
large as the peculiar velocity squared v 2 ~ (<9$) 2 . The amplitude of these terms are only 
the same size as $, with contributions from higher derivative terms being contained in 
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pure divergence terms which are zero by statistical homogeneity. Higher-order terms are 
suppressed by further factors of $. Any infra-red contribution must be a gauge effect, 
and can be renomalized away. And so, backreaction cannot be significant. 

Much of this argument is correct. Yet it also contains elements which are startling: 
second- order terms are the same size as first-order. In particular, the second-order 
potential has amplitude $( 2 ) ~ $, which is really not supposed to happen! Furthermore, 
if we calculate dimensionless measures of the time derivative of the Hubble rate we 
find contributions from perturbation theory like (<9 2( 3>) 2 , which are of order unity, and 
are actually divergent in the UV for a scale-invariant initial power spectrum; using 
any natural cutoff - from the end of inflation, say - gives colossal numbers. As we 
show here, these are not just a curiosity: we find the same thing if we calculate the 
monopole of the distance-redshift relation, a physical observable. We argue, as others 
have done [7,8,11,37,38], that it appears that higher-order perturbations will make 
things worse, not better. 

This is not to say that backreaction is necessarily large or significant, let alone 
responsible for an apparent acceleration effect. Rather we argue that until we have a 
better grasp of the consequences of these divergences, and higher-order perturbation 
theory in general, it is impossible to conclude that backreaction is small. Clearly, 
linear perturbation theory is correct in many respects. This success comes from 
analysing anisotropic and inhomogeneous perturbations, and observations of the CMB 
temperature anisotropy and the matter power spectrum, which all confirm it is basically 
correct. But backreaction concerns the homogeneous contribution of these perturbations 
at higher-order when they are calculated in terms of initial perturbations. Are they 
small? 

2. Measures of Backreaction 

In cosmology it is implicitly assumed that by statistical isotropy and homogeneity we 
simply construct a smooth FLRW 'spacetime' from an lumpy inhomogeneous spacetime 
by averaging over structure. This might, in principle, have a metric 

ds 2 ff = -dr 2 + 4 7 ydyMy J ' , (1) 

where r is the average cosmic time and a©(r) an averaged scale factor, the subscript D 
indicating that it has been obtained at a certain spatial scale V, which is large enough 
so that a homogeneity scale has been reached; in this case 7^ will be a metric of constant 
curvature. Unfortunately, we don't know how to construct ds 2 ff in an invariant way - 
there is no unambiguous averaging scheme. We don't know what field equations this 
metric would obey, nor do we know how to calculate observational relations; in fact we 
do not really know that any of the standard tools of general relativity survive. There is 
no reason to believe any of these would be as in GR, as it has only been tested in the 
solar system where the averaging problem does not apply. 



Is backreaction really small within concordance cosmology? 



4 



While different approaches to averaging exist which provide corrections to the field 
equations, those corrections require an underlying model with structure. In effect, it is 
as difficult to calculate a model with structure as the smooth average model which best 
corresponds to it. The scale factor itself comes from the solution to the field equations, 
but its first time derivative gives the Hubble expansion rate which is given algebraically 
from the Friedmann equation. Both of these change with averaging [1,3]. We can 
estimate what the averaged scale factor in Eq. (1) might be by different properties of it 
directly within the standard model, using perturbation theory. Before we get to that, 
let us first consider measures of the expansion rate locally in a spacetime; then consider 
averaging spatially; and finally let us consider how averages of observables look in a 
general spacetime - which is what we actually use to infer properties of a-p- 



2.1. Exact spacetime dynamics 

Given a 4- velocity u a , the local volume expansion is defined covariantly as 

6 = 3H = W a u a = 3u a W a \n£ = 3~ i , (2) 

where u a is the local fluid 4-velocity with length scale £. Locally, the expansion rate 
obeys the usual Friedmann and Raychaudhuri equations: 

Q^ = 3[p + A + a 2 -l 3 n] (3) 
= - § 9 2 - \{p + 3p) + A - 2a 2 + 2u 2 + A a A a + div A , (4) 

where the definitions of the terms are given in the Appendix. For an exact solution 
to the field equations, the Friedmann equation is a first integral of the Raychaudhuri 
equation, along with other equations such as the evolution equation for the shear - see 
the Appendix. From this we can define a local deceleration parameter 

q e = -1-3^, (5) 

= lL [UP + 3p) - A + 2a 2 - 2u 2 - div A - A a A a ] , 



e 2 

Up + 3p) - A + 2o 2 - 2u 2 - div A - A a A a 



(6) 



p + A + a 2 - 

which is commonly used as a dimensionless measure of the rate of change of the 
expansion rate. Over a small interval of cosmic time both and ge give important 
information about the dynamics of u a . 



2.2. Spatial averaging 

Let us assume for the moment an irrotational dust spacetime, and observers and 
coordinates at rest with respect to the dust. The average of a scalar quantity S may be 
(non-covariantly) defined as simply its integral over a region of a spatial hypersurface 
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V of constant proper time divided by the Riemannian volume (we follow [3]): 

(S(t, x)) v = ^- ( Vdet/id 3 ^ S(t,x) (7) 

Taking the time derivative of Eq. (7) yields the commutation relation 

d t (S) v - (d t S) v = (QS) V - (Q) v (S) v , (8) 

where is the local expansion of the dust, and we assume the domain is comoving with 
the dust. The dimensionless volume scale factor is defined as ax> oc V v 1/3 , which ensures 
(©)•£> = dtlnaD- Then, the equations of motion for this scale factor become: 

s(^)' = (p)v + ^-\lQv + (K) v } (9) 

3— + l( P )v = A+Qv, 

where Qt> = | [(0 2 )z> — — 2(cr 2 )x) is the kinematic backreaction term and 

°" 2 = \°~ab°~ ah is the magnitude of the shear tensor. The non-local variance of the 
local expansion rate can act in the same way as the cosmological constant, causing the 
average expansion rate to speed up, even if the local expansion rate is slowing down. For 
consistency, and assuming (p) v oc a^ 3 , these equations must lead to the integrability 
condition: 

-e-dtiQvat,) + —d t ((7Z} v a%) = . (10) 

A decceleration parameter q-p which describes the deceleration of the averaged 
hypersurface may be defined as: 

1 «£> / 1 -i \ 

<tD= ~ W — (11) 

l(p) I? -A + 2(a%-f[(0%-(0)| ) ] 



" (p) v + A + (o*) v - \ (K) v - 1 [{G^ - (0)!,] 

where — is given by the averaged Raychaudhuri equation. However, it should be noted 
that this is very different from the average of the local deceleration parameter, which is 

As we can see, it is not easy nor unique to estimate a© given a©, and we shall see that 
they are quite different in perturbation theory. 

These quantities come with a significant drawback: even for moderately sized 
domains, these quantities are unobservable because they are averaged on spatial 
hypersurfaces, whereas we only have observational access to our past lightcone (although 
see [39]). 
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2.3. Averaging physical observables 



Within the standard cosmology, cosmological parameters like the Hubble rate, and 
deceleration parameter are well defined based on the background metric. These 
cosmological parameters can be evaluated today, at t , by taking a Taylor series 
expansion of the scale factor a(t) 



a(t) = Oq 



1 + Hq (t _ tQ ) + I % Hi (t _ f )2 + 0([t- t ] 3 ) 



(14) 



Similarly, we can calculate the important observational quantities of interest, such as 
the relation between the angular-diameter distance and redshift d A (z): 



d A (z) 



z 

Ho 



l-±(3 + q )z + O(z) 2 



1 



z = H d A + - (3 + g ) (d A H o y + O (d A H y 



(15) 
(16) 



Of course, these series expansions are only valid in an FLRW background spacetime, 
but they can be generalised to an arbitrary spacetime using the Kristain and Sachs 
approach [40-42] . This gives a way to define a generalised Hubble rate and deceleration 
parameter in an arbitrary spacetime, which are observationally rather than dynamically 
motivated; that these have important differences was first shown in [43]. 
The corresponding equations to Eqs. (15) are: 



z = [K a K b V b u a ] d A + - [K a K b K c V a V h u c ] 4 + 



d A = 



1 



1 K a K b K c V a V b u c 



K a K b V a u b \ 
where the past pointing null vector is given by 

k a 



2 (K d K e V d u e y 



z + 



(u b k b ) 



-u a + e a , 



(17) 
(18) 

(19) 



which points in the direction e a at the observer, located at '0'. Using the notation 
whereby K An = K ai K a2 ■ ■ ■ K an and V ' An = V 0l V„ 2 • • • V a „, the n'th term K A "+ 1 V An u an 
for a given observer is an observable in principle, and will have a spherical harmonic 
expansion up to £ = n + 1. The monopole will be given by an all sky integral, 
J dVLK An+1 V An u an+1 , and is usually what we would compare to observations in an FLRW 
context. 

In an arbitrary spacetime, the first terms in these expansions are given by 

K a K b V a u b = |6 - A a e a + a ab e a e b (20) 
K a K b K c V a V b u c = l(p + 3p) + |0 2 - |A - lu a u a + a ab a ab + A a A a - fdiv A 

-e a [|D a 6 + |diva a + A a - |9A a - 1 {A b a ab 

+e {a e b) [E ab - \-n ab + \<da ab + u a u b + 1a c a a bc + e acd u c a b d - 2D a A b ] 
+e< a eV> [A a a bc - D a a bc ] (21) 
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The PSTF tensors e {Ae \ for t = 0, 1, 2, 3, in this expression are a covariant representation 
of the spherical harmonics when evaluated at a given point in spacetime (see the 
appendix). 

In the same way that we may evaluate various aspects of the spacetime by 
considering averages of the derivatives of a©, we can investigate how observations might 
be affected by the evolution of structure. For example, the observable K a K b V a Ub gives 
us a generalised Hubble rate which now has a dipole and quadrupole in addition to the 
familiar monopole term [43]. How shall we best compare higher order terms with the 
familiar FLRW expansion? Various definitions have been considered [6,43]. 

Note that the oIa(z) relation has powers of K a K b V a Ub on the denominator of each 
coefficient. Since this means dividing by a spherical harmonic expansion, it will be much 
simpler to work with the z(oIa) relation instead, and compare that with the FLRW 
relation. The definition we give below for the deceleration parameter is not unique, and 
depends on this choice. 

Let 

2 

#o° bs = E U ^ eAl = [K a K b V a u b ] (22) 
e=o 

where the %a ( are PSTF moments of the generalised Hubble rate H^ hs . Clearly the 
moments are simply, 

n = ^0, n a = A a and U ab = a ah , (23) 

evaluate at the observer. We now would like to relate the FLRW quantities (3 + qo)Ho 
to K a K b K C V V bU c , to define a generalised deceleration parameter. The simplest way 
is to define qQ hs via: 

(3 + g obs )^ 2 = [K a K b K c V a V b u c ] . (24) 

Then, writing qo bs = X^=o Q-Afi Al ■, where the PSTF tensors Qa ( are the multipole 
moments of the observational deceleration parameter, we have the monopole: 



Q= [|(P + 3p) - A + 6a 2 -2u 2 -2divA + 3A 2 } 



(25) 



Note that this is not the same as the local deceleration parameter associated with 6. 

We have, then, three different ways to analyse the second-derivative of the expansion 
rate, all of which have been written as variants of the familiar deceleration parameter. 
We have the local exact deceleration parameter of which we can calculate the average 
over a domain, the deceleration parameter associated with the deceleration of the 
averaged hypersurfaces, and the 'observed' deceleration parameter, which we have 
defined using the all-sky average (monopole) of the redshift-distance relation. In a 
general spacetime these are quite different things. 
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3. Backreaction in the standard model 

The standard model of cosmology ignores all the complexity of smoothing the spacetime 
by assuming that on 'large' scales (say larger than a few hundred Mpc) we can model the 
universe as homogeneous and isotropic. Linear and higher-order fluctuations describing 
structure then propagate as smooth fields on this background. In theory, when we sum 
over all these smooth fields we should end up with a metric which describes the universe 
we see, with structure properly described on all relevant scales (statistically speaking). 

We shall consider how backreaction comes about in the simplest cosmology which 
agrees with observations: a flat LCDM model with Gaussian scalar perturbations. In 
particular, we shall examine how backreaction affects the first and second derivatives 
of the scale factor, as given the Hubble rate and various deceleration parameters. 
Averaging FLRW perturbations has been discussed often in the literature: some 
authors investigate specifically the modification to the Hubble expansion rate or other 
variables [5,6,10-15,18,38,44-46]; others reformulate the average of the backreaction 
into an effective fluid [17,29,35,47-52], while one of the first attempts considered the 
important problem of how to calculate the averaged metric [53] . Many of these works 
consider only the case of backreaction in an Einstein-de Sitter model, in the hope of 
finding it responsible for dark energy. As it's more plausible that backreaction may lead 
to changes to the background at the level relevant for precision cosmology, we need to 
also investigate its effects in a LCDM model. 

3.1. Perturbation theory 

In the Poisson gauge to second-order in scalar perturbations the metric reads 

ds 2 = - [1 + 2$ + <5>^} dt 2 - aVid^dt + a 2 [(1 - 2$ - # (2) ) 7ij + h i3 ] dx l dx 3 , (26) 

The background evolution of the scale factor ait) at late times is determined by the 
Friedmann equation 

H{a) 2 ={^j = Hi [fi m a- 3 + 1 - Q m ] (27) 

where the Hubble constant H is the present day expansion rate, and Q m the normalised 
matter content today. The first-order scalar perturbations are given by \1/ (and are 
all that is required for observations at the moment), and the second-order by & 2 \ 
which are needed for a consistent analysis of backreaction. We also include the second- 
order vector modes Vi and tensor modes hij. In this gauge we have the metric in its 
Newtonian-like form, which we may think of as the local rest-frame of the gravitational 
field because it is the frame in which the magnetic part of the Weyl tensor vanishes 
when vectors and tensors are ignored [15]. 

For a single fluid with zero pressure and no anisotropic stress \1> = $, and $ obeys 
the 'master' equation 

$ + 4H<5> + A$ = . (28) 
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For a LCDM universe the solution in time to this is equation has <3> constant until A 
becomes important, and then starts to decay as A suppresses the growth of structure 
on all scales by about a factor of 2. We write it as <E>(i, x) = g(t)<&o(x) where g(t) is the 
growing solution to Eq. (28) normalised to g — 1 today (we can use g^ = g(t = 0) ~ 
|(3 + 2f2~ 045 ) as a very good approximation to its early time value). There is no scale 
dependence in the equation, which all comes from the initial conditions - usually a nearly 
scale-invariant Gaussian spectrum from frozen quantum fluctuations during inflation - 
and subsequent evolution during the radiation era. Evolution during the radiation era 
suppresses wavelengths which enter the Hubble radius compared to those which remain 
larger than it until the matter era begins. Consequently, in Fourier space, assuming scale 
invariant initial conditions from inflation, the power spectrum of $, V<s>, is independent 
of scale for modes larger than the equality scale, k eq = ^2Q m z eq HQ 0.07Q m h 2 Mpc -1 . 
A dimensionless transfer function describes the loss of power in the case of zero baryons 
(adapted from [54]): 

T (k) = In (2e + 0.134 «) (29) 

{ ' In (2e + 0.134 k) + [0.079 + ^f^] k 2 { ' 

where k = k/k eq . This is unity for k <C 1 and ~ (\xik)/k 2 for k 3> 1. The change in 
behaviour at the equality scale is important for backreaction because it is the modes 
larger than the equality scale which are primarily responsible for any backreaction at 
all. In essence, the equality scale determines the size of the backreaction effect. 
All first-order quantities can be derived from <3>; for example, 

v " = -sJhtA* + "*)' (30) 

is the first-order velocity perturbation, which governs the peculiar velocity between the 
matter flow and the rest-frame of the gravitational field. Meanwhile, the gauge-invariant 
density perturbation is 



Sp 2 

7 ~ 3H 2 n r , 



" 2 <9 2 $ - 3H ($ + H^j . (31) 



The second-order solutions for ^( 2 ) and $^ 2 ^ are given by [55]. These are 
complicated expressions involving time integrals over products of $ and its derivatives. 
For backreaction, however, the important thing is that if we take into account only 
terms which are important for backreaction we have simply 

= $(2) = B 3 (t)d- 2 d i d i (d i ^ d j ^ ) + B^d^od^o , (32) 

where B^{t) and B±(t) are time dependent functions given by integrals over g(t), with 
dimension Hq 2 . Their values today are given roughly by B 3 (t ) ~ 0.85/(f2^ 15 i^ 2 ) and 
B 4 (t ) ~ -0.25 /{Q^H 2 ). ^ and <5>^ also contain more complicated looking inverse 
Laplacian terms, but these drop out after ensemble averaging [15], so we ignore them 
here. We shall only explicitly calculate scalar quantities at second-order here, so we 
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don't require the exact form of the vectors and tensors; for reference, generically we 
have [56-59,59]: 

V^QdiQ, /iy (33) 
which we shall return to later. 

3.2. Backreaction of perturbations 

What is the backreaction of perturbations onto the the expansion rate and its associated 
deceleration parameter? Similarly, what are the corrections to the monopole of the 
distance-redshift relation from second-order perturbations? 

Expressions are quite cumbersome at second-order, so let us define the Hubble 
normalised dimensionless derivative 

d = {aH)~ 1 d. (34) 
At second-order, the expansion rate of the dust is 

i-(i + m-l(i + g)d 2 $ 

+ 1(3 - 2g) $ 2 - |$ (2) - \H~ 1 ^ + \ad k v\ 2) 

+ 27^T [I 1 + 2 9 + 9 2 ) + 9^ m (1 + g)] d k $d k $ (35) 

=>- H when evaluated today at the observer . 

This quantity is also the sky-averaged Hubble rate, measured by dust observers from 
observations on the past nullcone. We have used the fact that we are only interested in 
the growing mode solutions for $, and have written $ = g$ — gHQ where g = g/gH 
(a good approximation to its value today is g —1.31(1 — fi„ 4 ))- ^ n principle, then, 
the perturbative terms give the backreaction to the local expansion rate. Here, the 
divergence of the second-order velocity terms is 

an 2 m d k v k {2) = -^H ((2 - n m ) +4g + 2g 2 ) ($9 2 $ + d k <S>d k ^ (36) 

- 2 -n m (hB 2 ^ + S 2 ^^) + g) (<9 2 $<9 2 <£ + ~d k ~d 2 <b~d k §) . 

Note the terms with four derivatives in them - individually these are large, but they 
cancel out after ensemble averaging as we discuss below (or, in the language of [5], they 
only constitute a boundary term so should vanish on large scales). 

We may also calculate two measures of the deceleration within perturbation theory 
directly as they don't rely on averaging: ge and Q. These are rather messy, but have 
the same form up to coefficients. 

For example, the monopole of the deceleration parameter defined through the 
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distance-redshift relation becomes 
Q = -l + ^ m - 3(1 + £)(!- ft m )$ 



\9) 



9Q r 



(1 + 5) 



<9 2 $ 



-3g (4 + §ff) (1 - n m ) $ 2 - -^-(1 - Q m )(H^ + ¥V) 



2H 



1 



[3fi m (10 + 14$ + 6$ 2 ) - 8(2 + 5$ + 3$ 2 )] $<9 2 $ 
--^ [in^(13 + 12$) - 1^(2 - g - I2g 2 ) - ±(1 + 2g + g 2 )] ~d k ®~d k ® 



m 

4 



27fi 2 

m 



[l + 2g + g^didj^&Q 



(l + 2g + g 2 )d k $d'd k <S> 



27a 2 



27fi r 



[Q m (5 + 8g + 3g 2 ) - 4(1 + 2g + g 2 )] d 2 <S>d 2 <!> 



where it is understood that the rhs is evaluated at the observer, and is not a spacetime 
function in the usual sense. The other deceleration parameters are given in the appendix. 



3.3. Spatial Averaging 

Consider now the average of a variable over a spatial domain T>. The Riemannian 
average of a quantity T may be defined as, 



(T) 



1 



v 



Vdet/idVr 



v Jv 



(37) 



Here, the domain is on a hypersurface with 3- metric hij. If we choose that hypersurface 
to coincide with the spatial surfaces of the metric in the longitudinal gauge - the 
gravitational rest-frame - then this can be easily expanded in terms of the Euclidean 

average defined on the background space slices, (T) = / d 3 x T / / d 3 x, as: 



v 



v 



(T) grav = T(°> + (Y«> + (T< 2 >) + 3 [(TW)(*) - (T«*)] , 



(38) 



where T^ ), and T^ 2 ^ denote respectively the background, first order and second 
order parts of the scalar function T = + + Y^ 2 ). Note the important term 
in square brackets, which encapsulates the relativistic part of the averaging procedure. 
To link with Buchert's formulation above we would like to average in the hypersurfaces 
orthogonal to u a - the local rest frame of the dust - in which case we have instead 



(Y)x> = T(°) + (Y«) - <?/Y(°)($) + (Y^> 

+ (1 - 3H 9l ) [<Y«ct>> - (Y( 1 ))($)] - 9l {*fW) 

($ 2 ) + 3^(l-^)t (0) ($) 2 



(39) 



+91 



3HT m + -Y (0) 
2 



it(°) J\t' [($( 2 )) - ($ 2 > - 2<^%) - 2g I a- 2 (v\ 1) d&) 
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where gi = ^ f* g(t')dt'. So, for example, to find (q) v use = — 1 + §f2 m , 

TW = _3(1 + g)(l - n m ) $ + [(9 + 6g) - 4fi m (l + g)} <9 2 $ and T< 2 > = • • • all 
the other stuff • • • , in this expression (cf the Appendix). 
With this we now have 

(@)z> = the average of the expansion rate, averaged in the fluid rest frame 

(qe)v = the average of the deceleration parameter, averaged in the fluid rest frame 

q v = the deceleration of the scale factor associated with (0)x> 

Q = the monopole of the locally observed deceleration parameter (40) 

We now turn to estimating the amplitude of these quantities given standard initial 
conditions from inflation. 

3.4- Ensemble Averaging 

We define our Fourier transform as (suppressing any temporal quantities) 

$(s) = ^-L_ J d 3 fc$(fc)e*-, (41) 

where $*(fc) = $(— k). Note that we assume $(cc) is defined over an infinite volume 
for this definition to be valid, as it requires the crucial identity (2ir) 3 5(x) = J d 3 ke tkx . 
(This can instead be reformulated in a large box with periodic boundary conditions [60] . 
The periodicity of the boundary conditions is instead contained in 5-functions at k = 
here.) 

Inflationary models typically provide us with initial conditions for $ in terms of its 
correlations, which appear in the form of an ensemble average or expectation value: for 
example, the 2-point correlation function is just 



C(k, k') = $(fc)$(Jfe') , (42) 

where an over-bar denotes an ensemble average. If $ is Gaussian, then the distribution 
is given entirely by C{k,k'), and higher correlations are given in terms of C using 
Wick's theorem, and the ensemble average of odd numbers of $ is zero. The probability 
distribution function is just P[$(fc)$(fe / )] ~ exp ^c^k') • Statistical homogeneity of 
$(cc) implies that different modes are uncorrelated: C(k,k') oc S(k + k'); statistical 
isotropy implies that the proportionality function cannot depend on the direction 
of k [61]. Hence we have our power spectrum: 



2tt 



2 



$(fc)$(fc') = -^V^(k)5(k + fe') . (43) 
Assuming scale-invariant initial conditions from inflation, this is given by 



V*(t,k)=[^\ g(t) 2 T(k)\ (44) 
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where A^ is the primordial power of the curvature perturbation [62], with A^ m 
2.41 x 1(T 9 at a scale k CM B = 0.002Mpc _1 . 

What is the relation between ensemble and spatial averaging? Usually we assume 
ergodicity which roughly means here that, because the space is assumed infinite, the 
ensemble can be found within the space itself - i.e., any realisation exists somewhere. 
Consequently we should be able to replace ensemble averages with spatial ones, and vice 
versa. This only works if the spatial averages are taken over an infinite domain on the 
background, however, which is not necessarily what we are interested in when considering 
the averaging problem. We would like to know how properly averaged quantities change 
with scale. To see this, define the Euclidean average over a domain V using a window 
function W(x — x'). Then 

(X(x'))(x) = i J d 3 x'W(x - x')X(x') , (45) 

where V(x) = J d 3 x'W(x — x') gives the volume associated with W. For two correlated 
Gaussian random fields A and B satisfying A(k)B(k') = A(k)B(k)5(k + k') we have 
A(x)B(x) = (27r) -3 J d 3 kA(k)B(k) directly. However, when averaged over a finite 
domain, we have instead: 

(A(x')B(x'))(x) = ^ J d 3 x'W(x - x')A(x')B(x') 

= j^y 3 J d% J d 3 k 2 W(k 1 + k 2 ; x)A(k 1 )B(k 2 ) 

where W{k ; x) = V^x)- 1 f d 3 x'W(x — x')e tkx ' is the Fourier transform of W/V 
anchored at x. In the limit we average over all space, W/V — > 1 then W(k) = (27i) 3 5(k) 
and is independent of the anchor point x. Clearly we recover the ensemble averaged 
result in this case, provided the fields are isotropic. 



3. 5. The size of the backreaction 

To evaluate our various averaged quantities we could use a realisation of $ given an 
inflationary model. Alternatively, we can assume a spectrum for $ and evaluate the 
statistics of the quantity in question. This allows us to calculate the expectation value 
of averaged variables as well as their ensemble variance, in terms of integrals over the 
power spectrum of $ multiplied by powers of k. The reason we must go to second- 
order now becomes clear when we calculate the expectation value of an average: for 
Gaussian perturbations from inflation, the ensemble average of $ is zero, which implies - 
assuming ergodicity - that when averaged on the background over a very large (strictly, 
infinite) domain they are zero too. Thus, the second-order terms provide the principal 
backreaction effect; the first-order terms give the statistical variance. 

Let us estimate the approximate behaviour of each type of terms which appears. 
The relations for determining the scaling behaviour for the backreaction terms are {n+m 
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is even) 

_-Q(m+3n)/2 

(ojf) 

The inverse Laplacian term in <5>W satisfies (d^did^d^od^o)) = |(^$ 9^ ) [15]. 
We also have that, since $ is statistically homogeneous and isotropic [15] 



QmQ QnQ = ^_L__ J dkk^^V^k). (46) 



Q2<$>(2) = Q2^(2) = 0> Q kV k^ = o ? ( 47 ) 

which means that all the potentially large terms in the second-order Hubble rate don't 
contribute to the expectation value (see below). For the non-connected terms we have 

_ _ ( _])(m+3n)/2 poo 

= K > J dkk m+n - l W(kR D ) 2 V*(k), (48) 

where W is an appropriate window function specifying the domain. Typically this will 
become a delta-function as the domain tends to infinity (e.g., if the window function in 
real space is a Gaussian of width Rj>, then in Fourier space it is a Gaussian of width 
l/i?x>, centred at k — 0). Note that the connected terms have no dependence on the 
domain size or shape at all, and that the domain dependence arrises from the non- 
connected terms - these in turn come from using the Riemannian volume element. The 
integral can be written as 



— I d**—' P .(t)=(^) (£) / dnn m+n ~ 1 T(n) 2 . (49) 

Here, ku = Hq 1 is the wavenumber of the mode entering the Hubble rate today, and 
keq/ku = \f2VL m z eq ~ 40. Using z eq ~ 2.4 x 10 4 Q m h 2 implies the important relation 

2.4 Q 2 m h 2 . (50) 



k x 2 

k 



For pure CDM and a scale-invariant initial spectrum, the integral behaves as, replacing 
Jo°° ^ JTir V w here necessary: 

- ln(/tiR) for m + n = 
dKK m+n - 1 T(/t) 2 « { 3.9 for m + n = 2 (51) 

-^( K uv) f° r m + n = A 

The function F is roughly F(x) ~ 0.44 x 2AA for 1 < a: < 10, ~ 70or ai (log 10 x) 4 - 75 for 
x >> 1, and approaches ~ 53 In 3 a; as a; — )• oo. For integrals with the window function 
inside, W(K,k eq RT>), we can roughly replace Kuv ^ ^/Rvk eq , though this depends on 
the details of the window function used. Combining the above equations allows us to 
calculate reasonably precisely the size of each term at second-order. 

The first type of term, $ 2 , is nominally tiny, 0(1O~ 10 ). It also tells us that the IR 
divergence in $ 2 must be cut-off by hand, corresponding to the first mode to leave the 
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Hubble radius at the start of inflation. There has been speculation that this might lead 
to very important effects and even mimic dark energy [6,7,11], though this has been 
criticised [21-23]. Given that we only measure the primordial power spectrum to be 
nearly scale-invariant over a comparatively narrow range of scales, the appearance of 
($ 2 ) implies that it can't be too tilted to the red on super-Hubble scales, or else it would 
give a sizeable backreaction effect, though this might be unphysical. A red spectrum 
would convert the logarithmic divergence into a power-law one, and constraints on the 
largest mode would be much stronger (e.g., km/ha < 10~ 80 for n s = 0.95). 

The term primarily responsi ble for setting the fundamental amplitude of the 
backreaction in the Hubble rate is ($<9 2 $) oc (k eq /k H ) 2 - It is nominally quite small, 



(^ 2 $) 2 K q a 2 T t 



*k7TTr - ~ 10- (52) 



for the concordance model. (The overall effect is somewhat larger than this due to the 
contribution of several such terms.) This gives sub-percent changes to the Hubble rate 
from backreaction, though non-connected terms make it significantly larger on small 
scales [12]. Yet, we can observe that the backreaction is small because the equality scale 
is large in our universe, which is because the temperature of matter-radiation equality is 
very low. Modes which enter the Hubble radius during the radiation era are significantly 
damped compared to those which remain outside until after equality; so, the longer the 
radiation era, the less power there is on small scales to cause a significant backreaction 
effect. For a scale-invariant spectrum, then, it may be considered that the long-lived 
radiation era is the reason that the dynamical backreaction is small. The temperature 
will have to drop by several orders-of-magnitude before backreaction in the Hubble rate 
becomes significant. 

In the variance of the Hubble rate, the deceleration parameter q&, and the ensemble 
average of Q divergent terms appear: d 2 §d 2 &, which are of order the density fluctuation 
squaredf . The origin of such terms arises from taking the proper time derivative of the 
Hubble rate which contains spatial derivatives with respect to the coordinates we are 
using. The integral F overcomes the ~ 10~ 9 x (2Q m z eg ) 2 pre-factor around k uv ~ 10fc eg , 
so these terms are big, and difficult to know what to do with: here, the UV cutoff is really 
a measure of our ignorance. Within linear perturbation theory it should be set by the 
end of inflation and the reheating temperature, as well as the small-scale physics of dark 
matter, both of which are sub-pc scales today. Replacing the UV cutoff with a smoothing 
function in $ implies that we might do better in perturbation theory to smooth order- 
by-order, and calculate second-order terms from smoothed first-order ones, rather than 
calculating expectation values or spatial averages directly at a given order. Even for 

| While they appear in the Hubble rate through the second-order velocity perturbation, the ensemble 
average conspires to cancel them out for an initial spectrum which is homogeneous. Because the peculiar 
velocity appears as a pure divergence, the spatial average is a boundary term which is small on infinite 
domains due to assumed homogeneity - in our analysis this is contained in assuming ergodicity of the 
first-order $. 
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domains much larger than the non-linear scale, where linear perturbation theory breaks 
down (somewhere around a few Mpc), F is quite sizeable, and so we have backreaction 
terms 0(1) (note that F(l) = 0(1), and the pref actor is also 0(1)). From this, we also 
recover that the variance in the Hubble rate is 0(1) on scales of Mpcs. 

Observationally, can there be any signature of backreaction? When measuring the 
Hubble rate, perturbations are significant in the variance of the Hubble rate on sub- 
equality scales [12,14]. Perturbations affect the whole distance-redshift relation which 
has been calculated to first-order by [6,24,27]. Corrections to the luminosity distance 
include corrections <9 2 $, just like for the Hubble rate. This allows us to see that the 
variance includes divergent terms like (<9 2 $) 2 - the lensing term - which appear in 
the variance of the all-sky average of the luminosity distance - see [24]. We have 
discussed here how the problem at second-order is much more serious: the divergent 
terms appear in the expectation value of the distance-redshift relation itself, and not 
just in its variance. 

Consider the ensemble average of the monopole of gg bs . This should be the value a 
typical observer would expect to find - from nearby supernovae for example. Considering 
only the dominant terms with 4 derivatives in them, and using (k eq /k H ) 4 rs 5.7 fl^h 4 , 
we have (9(1) corrections to the observed deceleration parameter: 

q » -l + hi m + n 4 m h 4 [i.06 + 0.03(1 - n m ) 2 - 1.4(1 - n m ) u ] f(k vv ) (53) 

similarly , 

qe ~ -1 + + ^h 4 [-0.24(1 - fi m ) 3 - 94 + 0.665^f] F(«uv) (54) 

where the coefficient in square brackets are reasonable empirical estimates for fl m > 0.1 
(accurate to a % or so). In Buchert's definition of qx>, the divergence is nicely controlled 
by the domain size because terms (<9 2 $) 2 always appear as (d 2 <&y v . In that case the UV 
cutoff can be thought of as the smoothing scale. 

It is an important open question to find out what happens to the whole distance- 
redshift relation to second-order, where an overall ensemble averaged offset to the 
luminosity distance will be present. As we have shown, this will include terms (<9 2 $) 2 , 
which is certainly important enough to disturb the cosmic concordance. 

3. 6. Higher-order perturbation theory 

Would higher-order perturbation theory affect these conclusions? On the one hand 
it seems clear that higher-order perturbations should be suppressed. Provided <3> is 
Gaussian, only even orders will be important, once ensemble averages are taken. We 
might expect the largest terms at any order n to behave like <3>( n ) ~ (9$)™ (e.g., from 
relativistic corrections to the peculiar velocity), the ensemble average of which goes like 
A^(/c e(J //c#) n . Terms which appear in the Hubble rate at order n of the form 9 2 $^ do 
not have enough derivatives to overcome the suppression from (d&)( n ~ 2 ^ terms. By this 
argument, second-order should be as large as it gets, and backreaction from structure 
is irrelevant. 
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On the other hand, others [7,8,11,37,38] have argued that at higher order terms 
such as ((9$) 2 ((9 2< l ) ) n ~ 2 are the norm. Simply squaring terms like §d 2 Q gives problematic 
terms. In this case, from 4th order on perturbation theory needs a UV cutoff - at least 
as far as calculating averages is concerned. Even if not divergent, if (<9 2 $) n ~ 2 ~ 1 then 
higher-order terms are at least as large as at second-order and must be included to 
evaluate backreaction properly, and so correctly identify the background. But do these 
terms cancel out? 

Let us consider the expansion rate at fourth order, as an example. If we are 
interested in the expectation value, the leading terms, which contain 6 or more 
derivatives, are: 

tt(4) i o 

' " (1 + 3g + 3g 2 + g 3 ) d k $d k <f>d j d 2 $d j <f> (55) 



H 729m 



m 



[3(1 - g 2 - 2f) + tt m (27 - 67g + 30g 2 )] 4$<9 fc $<9 2 $<9 2 $ 



243fif; 



m 



1 fc. 

[2(1 + 3g + 3g 2 + f) + 7fi m (l + §)] d k ^d k d^d^d 2 ^ 



243^ 



_L_ [d k ~d 2 ^ - ~d k ~d 2 ^) {H~d k ^ + d k ¥ 2 ^ 



81^ 

+ terms of the form: $ (4) , d k v k {A) , d k <&d k §d 2 <& {2 \ 

<9 fe <9 2 $ (3) <9 fe $, d 2 §d i §d,M 2 \ d k $ {2) d k ¥ 2 \... 

Here, we have assumed the metric can be written in the longitudinal gauge up to 0(4), 
and we have ignored all stand-alone cubic terms as they do not appear in the ensemble 
average. We have also ignored coupling with vector and tensor modes which will give 
potentially important contributions. We have also grouped $( 2 ) and \I/( 2 ) together, 
schematically unless explicitly written out. Note that there is the potential for 8 
derivative terms from terms like (d k d 2 ^ 2 ^ — d k d 2 ^^)d k & ( - 2 ^ which appear. However, on 
closer inspection they cancel on small scales where $^ 2 ^ ~ ty( 2 \ so should not dominate 
the expression. We have included the coefficients where it is obvious how to calculate the 
ensemble average using Wick's theorem; where terms include $( 2 ) and $^ 3 \ it is more 
involved and will be considered elsewhere. The full calculation of is non-trivial, 
and we give this expression only to demonstrate what might happen. 

The amplitude of the 6 derivative terms may be estimated from Eq. (49). For 
example, we may estimate that 

d k $d k <bd 2 <&d 2 <& ~ d k <$>d k <& <9 2 $<9 2 $ ~ (^*) f^A x 3.9 x FU UV ) . (56) 

V 5#oc / \k H J 

Now, six powers of the equality scale overcome three factors of A^, using Eq. (50), 
making this at least the same size as first-order. If we take linear theory at face value, and 
assert that the UV cutoff should be on pc scales, then Kuv ~ 10 8 , and F ~ 10 5 ~ A^, 1 . 
Plugging in the numbers, we find these contributions to the Hubble rate - very naively - 
are of order unity. That is truly astonishing. This result is very sensitive to the UV 
cutoff, as well as the baryon fraction which helps reduce small-scale power. 
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If we consider what will happen to other quantities at fourth-order, it is clear that 
such six-derivative-with-four-Phi terms appear also from purely relativistic degrees of 
freedom. For example, in the Riemann tensor terms such as (d^hjk)) 2 will appear. Using 
the fact that ~ <&d(idj)<&, coupled with the fact that time derivatives give factors of 
k for gravitational waves, we see that such relativistic contributions can also be large; 
this cannot be associated with any Newtonian effect. Do they add up to something 
substantial? 

Perturbation theory suffers from higher-derivative terms because of the assumed 
scale-invariant primordial power spectrum, with mild damping of short wavelengths 
during the radiation era. Without this damping, even the $ 2 terms would lead to 
significant backreaction, as modes all the way up to the inflationary cutoff would be 
included. As it is, contributions for $ 2 and Qd 2 Q are cutoff by the Hubble scale 
at equality (k eq ), which is huge in comparison to more fundamental cutoffs. This 
suppresses modes shorter than the equality scale, but only in a power-law way, not 
exponentially. Yet it only takes 4 derivatives to overcome this power-law suppression at 
second-order: (<9 2 $) 2 . Exponential suppression of modes below around the Silk scale, or 
some hypothetical non-linear scale of about that size, would keep such terms just about 
under control [15]. 

One can argue that we should insert a non-linear cutoff scale k^L into our integrals 
over P$ to represent the effective reduction of power expected on small scales which 
hides the UV divergence. Indeed [35] claim a relativistic virial theorem provides such a 
scale. Of course, such a theorem cannot exist in general because energy is radiated to 
infinity in GR, and only the stationary part of a system virialises [63]. Nothing dynamic 
is truly isolated. This is realised in a non-trivial gravitational wave background from 
structure formation [58]. But this idea does imply that generically, at a fixed time, UV 
divergences represent something unphysical and can be renormalised away at that time 
(in [35] they deal with divergences by using a smoothing process on the background, 
which is effectively a UV cutoff). It is precisely this renormalisation which may lead to 
a significant backreaction. 

In integrals over the linear primordial power spectrum, however, seem to be the 
wrong place to insert a UV cutoff. While one might expect a reduction of power 
in the overall gravitational potential $( total ) today (if such a thing were meaningful) 
compared to $, with a change in behaviour at some scale k^L where linear theory 
becomes inaccurate, this should come about naturally from higher-order perturbation 
theory. If we actually knew the late time power spectrum everything would be fine - 
but an important aim is to calculate it directly from initial conditions, especially when 
we want to know how the background is renormalised by structure formation. In effect, 
we do not know what the relation of the background today is to the background at early 
times (if backreaction turns out to be an important effect). Higher-order perturbation 
theory presumably has to rely on the power at small scales in $ to calculate $( total ). 
There is no a priori reason to expect small scales to decouple from large order-by-order, 
even if they actually do decouple physically for 'virialised' regions after the series is 
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summed. So, an ad hoc insertion into integrals over $ doesn't seem to be the solution. 
It seems higher-order perturbation theory should be considered in detail to solve this 
problem. 

4. Discussion 

We have discussed here some examples of where the backreaction problem in the 
standard model of cosmology is non-trivial and not necessarily negligible. The two 
key examples we have discussed are the cases of the deceleration parameter which gives 
a dimensionless measure of the second time derivative of the scale factor, and the Hubble 
rate evaluated at fourth-order. In these cases perturbation theory does not give sensible 
answers because of terms which have significant power in the UV.§ 

In some respects dynamical backreaction is small by good fortune: the universe is so 
hot and had such a long radiation era that small-scale power is significantly reduced over 
its scale-invariant initial conditions (assuming those to be the case). The backreaction 
terms in the Hubble expansion, &d 2 &, are also the largest ones which appear in the lhs 
of the Einstein Field Equations, because the Einstein tensor has at most two derivatives 
of the metric in it. In some respects then this settles it: backreaction is small by virtue 
of there being a very small hierarchy of scales between the Hubble scale at equality and 
the Hubble scale today (they are only a factor of 50 apart in comoving terms). In this 
evaluation of backreaction, then, what happens on scales smaller than the equality scale 
is of little relevance. This is perhaps counter-intuitive given how we normally think of 
backreaction arising from small-scale structure - it is really power on very large scales 
which are responsible for the backreaction effect. 

On the other hand, however, we should be able describe the universe using a tetrad 
version of the Field Equations, or the covariant formulation given in the appendix. 
This is entirely equivalent to the EFE, but reformulates gravity as a system of first- 
order PDE's in the Ricci rotation coefficients and Weyl curvature tensor. In such a 
formulation higher derivative terms appear which are divergent implying backreaction 
may not be irrelevant. Furthermore, we should also be able to make sense of physical 
observables such as the fundamental distance-redshift relation, or the geodesic deviation 
equation. In these cases, terms such as (<9 2 $) 2 occur frequently, and as we have seen 
can appear in both their connected and non-connected forms giving rise to divergent 
backreaction terms, at least on small scales. 

In Buchert's interpretation of backreaction, where qx> describes the deceleration 
of the average scale factor, the divergence is neatly controlled by the domain size, 
and consequently has an elegant and straightforward interpretation, although it is not 

§ In [6, 64] it was claimed that a divergence in the IR can cause a large backreaction effect, 
whereas [21-23] have strongly disagreed with this idea. To make such terms large requires fluctuations 
out to enormous scales beyond the horizon. The UV divergence, by contrast, arises from the use of scale 
invariant initial conditions all the way down to small scales, which is usually found in many models of 
inflation. 
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observable. This also appears to be robust against various possible gauge effects such 
as averaging hypersurface [18]. The fact that there is a significant difference between 
the 'cosmological' deceleration and one which is smoothed on scales of a few Mpc, say, 
is expected since that is the scale where the Hubble flow kicks in. 

On the other hand, our two other definitions of the deceleration parameter - which 
do not depend on Riemannian averaging - reveal significant problems. Consider 
and qQ hs - the local deceleration parameter defined relative to the dust observers, 
and the observed one defined through the distance-redshift relation. If we calculate 
the expectation value of either of these we get enormous terms. Consider cutting 
off at a scale suggested by either a scale associated with the end of inflation, or 
from the dark matter suppression scale, which is around pc scales. Then we have 
k uv ~ lpc^/fdOOMpc)- 1 ^ io 8 . Our divergent integral then gives F(10 8 ) ~ 10 5 ! Of 
course, we have assumed a purely dark matter transfer function, and there is extra 
suppression from baryons on small scales, but this only reduces it by an order of 
magnitude for 20% baryons. One can reformulate the UV cut-off as a smoothing of 
the first-order potential to give something like the same interpretation one might give 
to Buchert's qx>, but this is rather ad hoc in this context - for (q@)x>, smoothing has 
already been done! Instead it maybe tells us that smoothing is necessary order by order 
in perturbation theory. That is, before constructing second-order perturbation theory, 
one necessarily must smooth structure below a certain scale. But why? Does this imply 
that the very notion of ergodicity needs to be made Riemannian: should ensemble 
averages be replaced by Riemannian spatial averages, rather than spatial averages on 
the background? 

This is very unsatisfactory. It is easy to devise universes where the magic number 
A^(/c eg //c_ff) 2 ~ 1 (though they may not be anthropically relevant). In an Einstein-de 
Sitter model, for example, one simply has to wait until such a time as H is small enough 
(though that doesn't happen with LCDM). If there were no radiation era things would 
be much worse. In such a case perturbation theory would simply break down at second- 
order and other formalisms for describing structure would have to be devised. The 
damping which happens at linear order to suppress a scale-invariant power spectrum 
would arise from non-linear structure formation, as fully non-linear GR is a well behaved 
theory. 

In our universe, where, apparently coincidentally, Aft(fc e9 /fc#) 2 ~ 1, it is 
A^keg/kn) 4 ~ 1 which seems to signal problems, because it is accompanied by a factor 
of F(kuv)- We can ignore this in many applications of perturbation theory because it 
simply does not appear. Should we be worried that they do appear in the examples 
calculated here? Perhaps not. Maybe there are ways around them, or perhaps they are 
some strange gauge artefact - though this seems difficult to realise since the distance- 
redshift relation is a physical observable (though its ensemble average is not). But we 
should be worried by the fact that (the ensemble average of) 2nd order perturbations 
are the same size as first, and 4th might be the same size as second - or larger if the 
(<9 2 $) 2 terms are to be taken at face value. It appears as if fourth-order perturbations 
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could give significant changes to the Hubble rate, even if one can dismiss the divergent 
examples we have discussed at second-order. Generically, we can see why (<9 2 $) 2 terms 
must appear at fourth-order because second-order potential terms are $<9 2 $; second- 
order gradient terms (such as v l , 2 C) are then ~ d§d 2 Q, and so have large expectation 
values when squared. It is tempting to assume that any terms with (<9 2 $) 2 in them 
all must cancel, or must be a gauge effect. However, as we have discussed, in models 
with <3><9 2( 3> ~ 1, we do not expect to have large terms cancelling just because they are 
large, so it need not be the case with (<9 2 <f>) 2 either. It is actually a combination of these 
terms together with scale invariant initial conditions which causes the problems. If they 
really do cancel in all possible models, at all orders in perturbation theory, this would 
be a very persuasive argument that backreaction is indeed small; if not, it is hard to say 
either way. 

Will perturbation theory even converge? Renormalization methods have been 
devised for Newtonian gravity [65,66] which is comforting, but only up to a point: vector 
modes which represent relativistic frame dragging have no Newtonian counterpart, are 
not that much smaller than $( 2 ) [59]; gravitational waves induced by linear scalars are 
larger and their time derivative is the same size as $( 2 ) [58]. These will easily mix 
together in any proper resummation scheme. However close Newtonian gravity may 
seem to relativity, GR is the theory of gravity we must work with to address these 
issues. 
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Appendix A. covariant formulation of the field equations 

The 1+3 covariant Ehlers-Ellis formalism provides a physically transparent formulation 
of the field equations full nonlinear generality (see [67,68] for reviews). The Ehlers-Ellis 
formalism is a covariant Lagrangian approach to gravitational dynamics, based on a 
decomposition relative to a chosen 4- velocity field u a . The fundamental tensors are 

h a b = 9ab + UaU b , S abc = 1] abcd U d , (A.l) 

where h a b projects into the instantaneous rest space of comoving observers, and e abc is 
the projection of the spacetime alternating tensor i] abcd = —^/—g^ia^b^c^d], and so 

r] abcd = 2u [a e b]cd - 2e ab[c u d] , e abc e def = 3\h [a d h b e h c] f . (A.2) 

The projected symmetric tracefree (PSTF) parts of vectors and rank-2 tensors are 

V( a ) = h a b V b , S( ab ) = ^h( a c h b ) d - -h cd h ab ^S cd . (A. 3) 
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The skew part of a projected rank-2 tensor is spatially dual to the projected vector, 
S a = \^abc^ bc \ and then any projected rank-2 tensor has the decomposition Sab — 
\Sh ab + e abc S c + S(ab), where S = S cd h cd . 

The covariant derivative V a defines 1+3 covariant time and spatial derivatives: 

= //' V,. /"-..,, I),./""..,, = ///'//", • • • h b f V d J e "... f . (A.4) 

The projected derivative D a defines a covariant PSTF divergence, divV = 
D a V a , div S a = B b S ab , and a covariant PSTF curl, 

curl V a = e abc D b V c , curl S ab = e cd{a WS b) d . (A.5) 

The relative motion of comoving observers is encoded in the PSTF kinematical 
quantities: the volume expansion rate, 4-acceleration, vorticity and shear, given 
respectively by 

6 = B a u a , A a = ii a , u a = curlw a , a ab = B {a u b) . (A. 6) 

Thus 

V b w a = ^Qh ab + e abc u c + a ab - A a u b . (A. 7) 

The PSTF dynamical quantities describe the sources of the gravitational field: 
the (total) energy density p = T ab u a u b , isotropic pressure p = ^h ab T ab , momentum 
density q a = —T^ a ' )b u b , and anisotropic stress n ab = T^, where T ab is the total energy- 
momentum tensor. The locally free gravitational field, i.e. the part of the spacetime 
curvature not directly determined locally by dynamical sources, is given by the Weyl 
tensor C abcd . This splits into the PSTF gravito-electric and gravito-magnetic fields 

E a b = C acbd u c u d , H ab = -e acd C cd be u e , (A. 8) 

which provide a covariant description of tidal forces and gravitational radiation. 
The Ricci and Bianchi identities, 

V [a V fe ]M c = R ab cdU d , V d C abcd = -V [a |i? 6 ] C - g%6]c}, (A.9) 

produce the fundamental evolution and constraint equations governing the covariant 
quantities. Einstein's equations are incorporated via the algebraic replacement of the 
Ricci tensor 

R ab = T ab _ _ T c g ab + Ag ab ^ (A.10) 

where T ab is the total energy-momentum tensor. 

The resulting equations, in fully nonlinear form and for a general source of the 
gravitational field, are: 
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Evolution: 




n + ( n + n )G + divo — — 2A a a — cr ab ir u 


(A 11) 


6 + l -Q 2 + hp + 3p) - A - div A = -a ab a ab + 2cu a u a + A a A a , 
o 2 


(A.12) 


Q(a) + ^©<?a + (p + p) A a + D a p + div 7r a = -a ab q b + e abc u b q c - A b ir ab , 
6 


(A.13) 


2 1 , 
W(«> + ^@w a + -curl A a = cr ab u , 

o 2 


(A. 14) 


2 1 

0~(ab) + 77@0"afc + E ab - -ll ab - D( a Aft) = -(7 c ( a (T fe ) C - UJ( a UJ b ) + A( a Aft) , 


(A.15) 


£ (a ft> + 9£ a6 - CMl\H ab +hp + p)(T ab + h(ab) + ^@7Ta6 + ^(a^) 

2 2 b 2 


(A.16) 


— —^-{aQb) + 2A c e c d( a H b ) d + 3o~ c ( a E b } c — u c e c d( a E b ) d — -o" c ( a 7r 6)c — 7^ c £' 


*r d 

cd(a^b) , 


f/^) + 9f/ a ft + curl E ab - -curl 7r af) = 3a c{a H b) - u e cd{a H b) 




3 1 c 
- 2A c e C(i ( a £ , ft) c - -uj( a q b ) + -o- c ( a e b ) cd q . 


(A.17) 


Constraint: 




div a; = A a w a , 


(A.18) 


div a a - curl w a - ^D a + g a = -2e: a ft c w fe A c , 


(A.19) 


curl (t„k A- D /,,ujh\ — fl,,h = — 2>1/«6l7m 


(A.20) 


1 l l 

div£ a + -div7r a - -D a p + -<dq a 

2 3 3 




1 3 

= £abcO~ dH c — 3H ab u + -o ab q — -e a 6cW q c , 


(A.21) 


div if a + ^curl g a - (p + p)w a 




= —£abc<J h dE cd — -e abc a b d 7T cd + 3E ab u b — -n ab u b . 


(A.22) 


The energy and momentum conservation equations are the evolution equations (A. 11) 


and (A. 13). The dynamical quantities p,p,q a ,n ab in the evolution 


and constraint 


equations (A.11)-(A.22) are the total quantities, with contributions from all dynamically 


significant particle species. 





Appendix A. 1. Moment Decomposition 

Here we give a summary of the covariant spherical harmonics we use, following [69]. 

An observer moving with 4— velocity u a at position x\ in a direction e a on the unit 
sphere (e a e a = 1, e a u a ) measures the luminosity of a distant supernova or galaxy. The 
direction e a can be given in terms of an orthonormal tetrad frame as, for example 



e a (9, 0) = (0, sin 9 sin 0, sin 6 cos (f>, cos 9) 



(A.23) 
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There is a 1 — 1 mapping between all symmetric trace-free tensors of rank / and the 
spherical harmonics of order £ [69-71]. Consider a spherical harmonic decomposition of 

oo oo +£ 

/ = E/' = EE Fp(x*)Y e m (6, 0) (A.24) 
1=0 1=0 m =-e 

where Y£ m (#,0) are the usual surface spherical harmonic. An alternative way to carry 
out this harmonic expansion is to write / in the form 

oo 

f = J2 Fa ^ Ai =F + F ^ + + F abce a e b e c + F abcd e a e b e c e d + .... (A.25) 

£=0 

where the spherical harmonic coefficients F A( are symmetric, trace-free tensors 
orthogonal to u a : 

F Ai = F {Ae) , F Aeab h ab = 0, F Aea u a . (A.26) 

Round brackets "(••)" denote the symmetric part of a set of indices, angle brackets "(••)" 
the PSTF part of the indices. We use the shorthand notation using the compound index 
A £ = ai a 2 ...a £ , and e Ae = e ai • • • e ai . The PSTF part of e<^> [70-73] is 

[1/2] 

e (A e ) = B ek h {aia2 ■ ■ ■ h a2k - ia2k e a2k+1 e a2k + 2 ■ ■ ■ e ae) , (A.27) 

fc=0 

where B ek are given by [70], 



{-l) n £\(2£-2n-l)\\ 
en ~ (£-2n)\(2£-l)\\(2n)W ( 8) 



Here [£/2\ means the largest integer part less than or equal to £/2, and £\\ = 
2)(£ — 4) • • • (2 or 1). The normalization for e Ae is given by [73], for odd and even £ 
respectively: 

— I e A2 ^oin = 0, and — / e A2e dfl = — *— , (A.29) 
AnJ 4w 4irJ 4w 2£ + 1 1 } 

which implies 

/ e Ae e Bm dn = — h {AlB ™\ (A.30) 

£ + m + l 

if £ + m is even, and is zero otherwise. The orthogonality condition for e^ A ^ then follows: 
J dQe^e {Bm) = 5 l m A e h^ {bl ■ ■ ■ h^ be) with A, = ^ + ^ ±^ , (A.31) 
which then implies the relation between / and its spherical harmonic moments f Ae : 

F Ae (x i ) = Aj 1 [ dQe {Ae) f(x\e a ). (A.32) 

Jilt 

Note the recursion relation 

f 2 

e (A e+1 ) _ e (a e+le (A e )) ^ /^K+i^g^-i}) (A. 33) 

relates the (£ + 1) — th term to the £ — th term and the {£—!)— th term. 
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We give here the full deceleration parameters for the cases we consider. First the 
deceleration parameter associated with the fluid expansion rate 6: 



q e = -1 + -tt m - 3(1 + g)(l- On)$ + 



9a 



"(I +9) 



<9 2 $ 



-33 (4 + \g) (i - n m ) $ 2 - — (i - n rn )(H^ + vi/( 2 )) 
+^(i - 3n m ) (H~d k v\ 2) - d k v k (2) ) - 1 (W 2 ) - dV>) 



1 



2701 



[3Q m (U + Ug + Qg 2 ) - (16 + 3Qg + 2g 2 )] $«9 2 $ 
[f n 2 m (ll + 12g) - 24fi(l - ^) - 4(1 + 2£ + g 2 )] d k <f>d k <f> 



27W 



+ 



(1 + 2g + 3 2 )<9 fe $<9 fe <9 2 $ 

[n m (5 + 8g + 3£ 2 ) - |(1 + 2£ + £ 2 )] d 2 $d 2 <5> . 



27Vl m 

The deceleration parameter an observer would measure from the all-sky average of the 
redshift-distance relation: 



Q=-l + -n m - 3(1 +g)(l- On)$ + 



\9) 



9a 



[1 + 9) 



«9 2 $ 



-3g (4 + lg) (1 - n m ) $ 2 - ^(1 - O^tf^ 2 ) + ^) 



1 



3H 



[3fi m (10 + 14^ + 6£ 2 ) - 8(2 + 5g + 3g 2 )} $<9 2 $ 



[in^(i3 + 123) - |On(2 - 5 - vg 2 ) - U l + 2 9 + f)] 

A ____ & 



2701 



27a 



(1 + 2g + g^did^W® - ^-(1 + 2^ + ^)d fc $c> 2 c> fe $ 
[ft m (5 + 8^ + 33 2 ) - 4(1 + 2g + # 2 )] «9 2 $<9 2 $ 



Finally, we give the deceleration parameter defined via Eq. (11), averaged in the local 
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rest-frame of the dust observers: 



qd = -i + -n m - 3 [(1 + g)(l - Q m ) + p m giH] ($) + 



2 $> 



+ {3 (4 + g - \?) + 9 9I H [(1 + g)(l - fim) + - ^ (f - fim) } 2 ) 

-3 {(4 + 5g) + |Q m (3 + 4£ + £ 2 ) - 3 9l H [(1 + $) (l - \Q? m ) - Q m (l + 2g + A 9l H)] 

-f^^ 2 (l-|^)}W 2 



f} 2 



[§f&(3 + 4$) - |Q m (2 + 5 - £ 2 ) - A(i + 2 ^ + £ 2 )] (cW$) 



+ -|(17 + 20^ + 6^ 2 )- 

+ 9I H (1-2$)- 
3 



3a 



l0 + 14^ + 4^ 2 ) 



+ 



2H 

9 

8 



3r2 m 

I -fi m )(tf($ (2) > 



l + 2^)-fi m (2 + ^) 



($)0 2 $) 



+ |a(2-3fi m )0 fe ^) + K0 2$(2) ) 

y'd*' ($ (2) ) - K$ 2 } - - 9j ^-{v\~d^) 

[f (1 + 2$ + g 2 ) - Q m (5 + 8g + 3g 3 )} 2 $) 2 . 



27Q 

Note that there are no connected (<9 2 $) 2 terms in this expression. 



(B.l) 
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